Detecting a physical difference between the CDM halos in simulation and in nature 
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Numerical simulation is an important tool to help us understand the process of structure forma- 
tion in the universe. However many simulation results of cold dark matter (CDM) halos on small 
scale are inconsistent with observations: the central density profile is too cuspy and there are too 
many substructures. Here we point out that these two problems may be connected with a hitherto 
unrecognized bias in the simulation halos. Although CDM halos in nature and in simulation are 
both virialized systems of collisionless CDM particles, gravitational encounter cannot be neglected 
in the simulation halos because they contain much less particles. We demonstrate this by two nu- 
merical experiments, showing that there is a difference on the microcosmic scale between the natural 
and simulation halos. The simulation halo is more akin to globular clusters where gravitational en- 
counter is known to lead to such drastic phenomena as core collapse. And such artificial core collapse 
process appears to link the two problems together in the bottom-up scenario of structure formation 
in the ACDM universe. The discovery of this bias also has implications on the applicability of the 
Jeans Theorem in Galactic Dynamics. 
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INTRODUCTION. 

Thanks to rapid development in numerical techniques, 
the N-body simulation method has scored considerable 
success in furthering our understanding of the formation 
and evolution of the large-scale structures in the uni- 
verse. But the simulations based on the CDM models 
has not worked so well on small scales, one aspect is 
the " cusp problem" : The central density profile of the 
simulation halos (Fukushige & Makino 1997; Moore et 
al.1999; Ghigna et al. 2000; Navarro et al. 1996) are 
too cuspy, and are inconsistent with the observations 
on dwarf galaxies and low surface brightness galaxies 
(Gentile et al. 2004; Primack 2003). It is also incon- 
sistent with the distribution of dark matter in galaxy 
clusters found from X-ray and gravitational lensing stud- 
ies (Katayama et al. 2004; Tyson et al. 1998 ). An- 
other aspect is that 10 ~ 100 times more substructures 
are predicted than is observed (Willman et al. 2002; 
Bilker et al.2005); — the so-called "substructure prob- 
lem". Some alternative models were introduced to deal 
with these problems, including the self-interaction dark 
mattcr(SIDM) model, the warm dark matter model etc. 
(Spergel & Steinhardt 2000; Gotz et al. 2003; Boehm et 
al. 2002). 

Bere we realize that our understanding of these two 
problems can be greatly helped by studying a common, 
basic technical difficulty: the huge difference between the 
number of particles used in simulations of CDM halo, 
(Nhaio) and the number occurring in nature: even in 
the most powerful "high-resolution" simulations today, a 
large dark matter halo contains only about 10 6 -10 7 par- 
ticles. Since we have to keep to a fixed mean density of 
the universe, the mass of each particle in the simulation 
has to be some 70 powers of ten the GeV candidates in 



particle physics (Bertone et al. 2001; Bertone et al. 2005) 
It is hard to accept that the mass of the CDM particle 
can be as high as the mass of a small galaxy. We have 
unavoidably introduced a physical bias in our simulations 
which would impact on small scale problems. The bias is 
that we have simultaneously reduced the particle number 
density n and raised the mass of each particle m by a huge 
factor, /ae ~ 10 70 . CDM particle systems with such a 
physical bias can be expected to have very much biased 
dynamical and statistical properties, and the clustering 
and evolution scenario of CDM halos we learn from them 
might be different from reality. This problem should be 
discussed in more detail. 



RELAXATION AND DETECTING METHOD 

Bow will such bias affect the simulation results ? One 
effect is it changes the mean free path of two-body re- 
laxation, L s . For a virialized CDM halo with mass M 
and radius R, one can roughly estimate L s as (Xiao et 



al. 2004): 
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where ~p = M/(4irR 3 /3) is the mean density of the halo, 
and p is its local density. For a virialized halo in nature 
that is composed of a huge number of GeV particles, it's 
easy to see from eq.(l) that L s R and that we can com- 
pletely disregard the two-body relaxation. Bowever, for a 
virialized halo in simulation, composed of a much smaller 
number of huge mass particles, we have L s ~ R, and the 
corresponding two-body scattering cross section is close 
to the value expected in the SIDM model. So, in the 
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FIG. 1: Time evolution of the distribution function f(rjm,t) and f(rjLi,t) of a stable and spherical halo. Different from the 
expected delta function of a halo in nature, the f(rj,t) of the simulation halo in 0.001 Gyr(dash/green line), 0.01 Gyr(dot/red 
line) and 0.1 Gyr(solid/black line) is becoming more and more dispersive. 



simulation, the effect of relaxation can no longer be ne- 
glected. Considering that star clusters that suffer similar 
relaxation can have their dynamical and statistical prop- 
erties seriously affected, leading to core collapse, evapo- 
ration, etc. (Binney & Trcmainc,1987), the discovery of 
this physical difference between CDM halos in simulation 
and in nature can help us understand the two small scale 
problems encountered in our current simulations. 

In an effort to quantify and detect the relaxation pro- 
cess we have designed two numerical experiments for a 
spherical CDM halo. Different from the traditional sta- 
tistical method (Binney & Knebe 2002; Dicmand et al. 
2004), we study the relaxation process from the view- 
point of the individual particle, and we find that the two 
integrals of motion, the energy and angular momentum of 
the ith particle of the halo, Ej,(x, v) and Lj(x,v) = |L|, 
can be effective detectors of the relaxation effect caused 
by the physical bias. 

Although, a globular cluster and a stable spherical 
CDM halo in nature are both composed of collision- 
less gravitational interacting particles, their values of 
L s /R ~ N show that their gravitational potentials on 
small scale are very different. For a CDM halo in na- 
ture, gravitational encounter can be neglected and the 
halo potential <p(r) is stable and spherical for each CDM 
particle. Then Ei and Li of each particle, being both 
integrals of motion, will be constant in such a system. In 
contrast, for a globular cluster or a CDM halo in simula- 
tion, L s ~ R, and this means the collisionlcss particles in 
these systems have to experience gravitational encoun- 
ters, and the potential <f> is neither spherical nor stable. 
Then the Ei and Li of each particle are no longer inte- 



grals of motion and are soon modified by encounters. — 
(The energy E and angular momentum L of the whole 
system are still the same. This is an important require- 
ment for the reliability of many simulation results, but 
it is not relevant to the problem on hand.) From this 
point of view, the CDM halo in simulation is more like a 
globular cluster than a CDM halo in nature. 

As measures of the relaxation effect in simulation CDM 
halos, we now define two parameters, 
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where AJSj = Ei(t) — E^ and ALi = Li(t) — £;o ar e the 
variations of Ei and Li from their initial values Eio and 
LiQ. For a CDM halo in nature, the distribution functions 
will simply be the time-independent delta-functions; for 
a stable, spherical halo in simulation, in contrast, they 
will be very different and will evolve in time. 



SIMULATION AND RESULT OF NUMERICAL 
EXPERIMENTS 

Here we present two controlled numerical experiments 
using the publicly available N-body tree code, GADGET 
(Springel et al., 2001), to detect the relaxation effect 
discussed above. Following Ma et al. (Gozdziewski et 
al.2005; Ma & Michael 2004; Michael & Ma 2004), we 
generate a stable spherical virialized CDM halo with a 
given NFW density profile (Navarro et al.1996): 



p{r) = p crit 6/[x(l + xf 
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FIG. 2: Time evolution of Ni e f t /N a ii for Ei (solid/black line) 
and Li (dot /red line). We can see the Ei and Li have been 
changed for most of the particles in a short time, different 
from the expected dash line of the halo in nature that will 
not be affected by such relaxation. 



FIG. 3: Time evolution of the distribution function f(j]u , t) of 
a spherical evolving halo. Similar with the stable halo in Fig 
1, the f(r),t) of the simulation halo in 0.001 Gyr(dash/green 
line), 0.01 Gyr(dot/red line) and 0.1 Gyr(solid/black line) is 
becoming more and more dispersive. 



where x = r/r a ,S = 200c 3 /3[Zn(l + c) - c/(l + c)], and 
c = r V i r /r s is the ratio of the halo's virial radius to scale 
radius. In all the runs, we set r s = I6kpc/h, the total 
mass of the halo M = 10 12 M@, and we use a force soft- 
ening of e = 2kpc/h. 1 We drew the particle velocities 
from a local isotropic Maxwellian distribution with the 
radius-dependent velocity dispersion computed from the 
Jeans equation. 

In the first experiment, we start with a stable halo of 
2 x 10 4 particles. We first let the individual particles 
of the NFW halo evolve for 5.0 Gyr, and we use the 
evolved halo as our initial condition, noting down the 
values Eio and of each particle at this time. Then 
we record the evolution of Ei and Li of each particle and 
obtain the distribution functions f{r)Eh t) and /(r?z,i, t) at 
chosen epochs, t. Theoretically, as integrals of motion, 
the Ei and Li of each particle of such a halo in nature 
will be constant and the distribution function f(rj) should 
remain a delta function all the time. However, as we can 
see in Fig 1, both /(^lJ and f{r]Ei) are spreading out as 
time goes on. The spreading means that within a period 
of 0.1 Gyr, more and more particles have experienced 
changes in their Ei and Li. 

We define Ni e f t as the number of particles which satisfy 
\f]Ei\ < 0.01 or \r]Li\ < 0.01 as measures of the evolution 



1 Our test runs showed that different values of p(r) and e etc. do 
have a small effect on the results, but will not change the general 
qualitative picture. 



of the distribution functions. As we can see in Fig. 2, 
for most particles in the simulated halo, Ei and Li were 
changed very soon, and that the changes were faster in Li 
than in Ei. In contrast, the value of Ni e ft/N a u remains 
constant at 1.0 for a CDM halo in nature that consist of 
a huge number of GeV particles. 

In the second experiment, we consider a spherically 
evolving halo. We multiply the velocity of each particle 
in the halo of the first experiment by a factor of 0.9 as 
the initial velocity. Then the halo will still be spherically 
symmetric, but it is no longer stable. In this case, Ei will 
be changed by the violent relaxation. Meanwhile, for a 
CDM halo in nature including numerous GeV collision- 
less particles, the Li of each particle will remain as an 
integral of motion and will not change. Fig 3 shows the 
distribution function f(fjLi) of the simulation halo. Simi- 
lar to Fig 1, the Li of each particle changed rapidly due to 
this hitherto unsuspected process of galactic relaxation. 

Finally, we tried to assess the effect of particle number 
on the distribution functions. We generated three sys- 
tems with, respectively, 2 x 10 4 , 2 x 10 5 and 2 x 10 6 par- 
ticles, the other parameters remaining the same as for the 
stable halo of the first experiment. Fig 4 shows the time 
evolution of ai h = (r^). We find that in all cases a^ L 
increases as a power-law function of time (within about 
O.lGyr), and, at any given time, a^ L is a fast decreasing 
function of N. In other words, the relaxation effect, or 
the bias, is greatly reduced when N is increased. 
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FIG. 4: a^ L = (r]^) as a function of time, for stable spher- 
ical halos with different particle numbers: 2 x 10 4 (dash/red 
line), 2 x 10 5 (solid/green line), 2 x 10 6 (dot/black line) and 
the spherical evolving halo with 2 x 10 (dash dot /blue line) 
particles. 



DISCUSSION AND CONCLUSION 

What is the difference ? — Our demonstration with the 
two numerical experiments brings out a physical differ- 
ence in microcosm between dark matter halos in simula- 
tions and in nature: although they are both composed of 
collisionless CDM particles, the process of gravitational 
encounter can be neglected in natural halos, but not in 
simulation halos. 

How will it affect the simulation results? — On small 
scale, the relaxation can greatly affect the behavior of 
individual particles in simulation halos. Take the case of 
the stable halo as an example. For the halo as a whole, 
the total energy E and the total angular momentum L are 
each conserved. What's more, statistically <fi(r) is still a 
static and spherically symmetric potential in a short pe- 
riod. However Fig 1 and Fig 2 clearly show changes of 
Ei and Li caused by gravitational encounters in a micro- 
statistical sense. In other words, from the point of view 
of the individual particle of the simulation halo, the po- 
tential is nolonger spherically symmetric or static. To be 
exact, for a halo that comprises fewer than approximately 
10 6 particles, gravitational collision must be taken into 
count. The halo is no longer a system of collisionless 
particles, and the collisionless Boltzmann Equation is no 
longer applicable. 

Will the general dynamical behavior of the whole sys- 
tem in the large scale statistics respect, be affected by 
such microcosmic effect? The answer is obviously yes. 
For the CDM halo in simulation, the Ei and Li of each 



particle can vary, which means the distribution function 
f(Ei,Li) of the whole system can be different from the 
natural case. And if we can't trust the f(Ei, Li) in sim- 
ulation, then the density profile p(r) and velocity dis- 
persion of the halo we learn from the simulation will be 
subject to doubt. Actually, whether we can express the 
distribution function /(x, v) in the form of f(E, L) needs 
to be considered (see the following discussion about the 
Jeans Theorem). There isn't any mature theory dealing 
with the effect of such microcosmic relaxation on p{r) so 
far, and more research is expected. 

However, by comparing our results above with a previ- 
ous research on dynamic behavior of stellar clusters (Bin- 
ney & Trcmaine 1987), we can learn a lot about how the 
difference caused by the particle number will affect the 
simulation results. Much research has been done on stel- 
lar systems (such as globular clusters) which comprise 
similar numbers of collisionless particles under pure grav- 
itational interaction. The dynamical effects caused by 
the process of gravitational encounter in stellar systems 
can also appear in simulation halos. A good example is 
the core collapse process: due to star-star encounters in 
the core, some of the stars get very high energies and 
evaporate away, while the core shrinks as required by en- 
ergy conservation. Following the way of studying stellar 
systems, we can also estimate the time to core collapse 
r of one CDM halo (Binney & Trcmaine 1987, eq. 8-72) 
as: 
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Where M and are the total mass and median radius 
of the halo, m is the mass of each particle of it, and t r h is 
the median relaxation time of the system . Now we can 
see, the core collapse time scale of an isolated spherical 
simulation CDM halo with M ~ 10 12 M Q ,r h ~ lOfcpc and 
including N ~ 10 4 particles (m = M/N ~ 1O 8 M ), will 
be about the age of the universe. Based on the estimation 
above, one might expect the CDM halo in simulation can 
avoid such core collapse process. Unfortunately, consid- 
ering: 

(1) According to the hierarchical structure formation 
scenario in a ACDM universe, small structures come into 
being first and large halos comes from the merging of 
small ones. 

(2) In the same simulation, with the same value of m 
above but much smaller M and r/j, eq.(4) indicates small 
halos will have a much shorter core collapse time scale 
than the age of the universe(one can also see from Fig. 4 
and eq.(l)). So these small halos in simulation will have 
to suffer the core collapse process, and form one artificial 
cusp center. On the contrary, such core collapse will not 
happen in a CDM halo in nature which includes very 
large number of small m ~ GeV particles, for they have 
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a time scale r much longer than 10 10 years. 

(3)Numerical studies (Hayashi et al.2003, Michael & 
Ma 2004) have shown the cuspy centers can survive from 
merging process of both the major mergers of galaxy ha- 
los and the evolution of substructures. 

Now we can figure out how do the bias affect the sim- 
ulation results: In a ACDM universe, small halos come 
into being first and the unexpected core collapse process 
lead to the artificial cuspy centers of them. Then these 
cuspy centers survive from the merging process later and 
contribute to both the "cusp problem" and the "sub- 
structure problem" . No matter if there is other explana- 
tion for these problems, the contribution of this effects 
must be taken into account when we analyze simulation 
results on small scales. One can even predict, with the 
improvement of simulation technology, the " substructure 
problem" can be even more serious if the particle number 
increased by only one or two orders. 

Last but not least, our discovery of the bias directly 
tests the premise of the commonly used Jeans Theorem 
in galactic dynamics. According to the Jeans Theorem, 
in a static and spherical symmetric potential well, we can 
express the distribution function /(x, v) in phase space 
as f(E, L) only when the and Li of each particle are 
all integrals of motion. Our result shows that, for the 
collisionless particle system with pure gravitational in- 
teraction which has a particle number N ^ I0 6 , gravi- 
tational encountering can not be dismissed. As a result, 
the Ei and Li for most of the particles will be changed in 
a short period of time (see Fig. 2) . For such systems like 
CDM halos in simulation and globular clusters, the Ei 
and Li are no longer integrals of motion and the Jeans 
Theorem is no longer applicable, while it is still valid for 
CDM halos in nature that are made up of a huge number 
of GeV particles. 
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